Diffusion dynamics on multiplex networks 
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We study the time scales associated to diffusion processes that take place on multiplex networks 
(a set of networks structured in interconnected layers) . To this end we propose the construction of 
a supra-Laplacian matrix, which consists of a dimensional lifting of the Laplacian matrix of each 
layer of the multiplex network. We use perturbative analysis to reveal analytically the structure 
of eigenvectors and eigenvalues of the complete network in terms of the spectral properties of the 
individual layers. The spectrum of the supra-Laplacian allows us to understand the physics of 
diffusion-like processes on top of multiplex networks. 
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Modern theory of complex networks is facing new chal- 
lenges driven by the necessity to understand properly the 
dynamical evolution of some physical processes on real 
systems. Among the open problems that are focussing 
the interest of the scientific community one of the most 
representative concerns the topological and dynamical 
characterization of what has been called multiplex net- 
works. The standard approach to the problem assumes 
a simplex view of networks where every edge (link) is of 
the same type and consequently considered at the same 
temporal and topological scale P. This is clearly an 
abstraction of any real topological structure and repre- 
sents either instantaneous or aggregated interactions over 
a certain time window. Therefore, to understand the in- 
tricate variability of real complex systems, where many 
different time scales and structural patterns coexist we 
need a new scenario, a new level of description Q. 

A natural extension which allows to overcome previ- 
ous drawbacks is to describe a network as a set of cou- 
pled layered networks (multiplex network) where each 
layer could have very particular features different from 
the rest, and in this way, define a richer structure of in- 
teractions. Some works have already focused on the de- 
scription and analysis of interconnected networks 
however mathematically grounded results about general 
dynamical processes running on them are yet to come. 

In this Letter we focus on a particular setup of multi- 
level networks in which nodes are conserved through the 
different layers of the multiplex (see Fig. [IJ . The current 
study analyzes a diffusion process that takes place at the 
whole system level, i.e. within and across layers. This 
setup could account, for instance, for a diffusion dynam- 
ics taking place on top of a social network of contacts. 
Admittedly, the latter is a network of networks, i.e. the 
aggregate of many different social circles or subnetworks, 
each having its own temporal or structural patterns (for 



example, think of our online activity which includes dif- 
ferent social networking sites such as Facebook, Twit- 
ter, etc.). The same applies to multimodal transporta- 
tion networks Q , on top of which individuals "diffuse" 
within and between different layers (e.g. bus, subway, 
etc.). Let us remark, however, that our interest here 
is not to solve a specific real problem but to illustrate 
the analysis of diffusion processes on top of these struc- 
tures. We propose a mathematical setting that allows 
to scrutinize the emergent diffusion time scales in multi- 
plex networks. We concentrate on diffusive processes, as 
they constitute a good approximation for different types 
of dynamical processes (e.g., synchronization and other 
nonlinear processes amenable of linearization j9[) whose 




FIG. 1: Example of multiplex network with M = 2 layers. 
Nodes are the same in both layers. The connectivity at each 
layer is independent of each other, the connectivity inter-layer 
is from each node to itself (dashed links). 
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dynamical properties can be captured by the behavior 
of the eigenvalues of the Laplacian matrix. For instance, 
the time needed to synchronize phase oscillators in a net- 
work is related to the second smallest eigenvalue of the 
Laplacian, A2 , and the stability of the synchronized 
state is determined by the eigenratio A7V/A2 [ilj . The 
spectral analysis of complex networks constitutes then 
still a promising area of research J1 2 L | l 3l|. Following a 
perturbative analysis of the spectra [1J|, our results al- 
low to get new physical insight about diffusion processes 
through the analytical determination of the asymptotic 
behavior of the eigenvalues of the Laplacian of the mul- 
tiplex (supra-Laplacian) when the diffusive coupling be- 
tween layers is either small or large. 

Let us consider a set up in which the diffusive dynamics 
is linearly coupled within nodes in each layer K, through 
a diffusion constant Dk, and among nodes in different 
layers K and L, in this case with a diffusion constant 
Dkl- The network at each layer is assumed to be con- 
nected and undirected, but it can be weighted. The state 
of each of the N nodes is represented as a vector indexed 
by layers xf (t) where the subscript stands for the node 
and the superscript for the layer. The equations describ- 
ing the dynamical evolution of the states of the nodes, 
considering a multiplex of the M layers, are: 



The dynamical properties of the system can then be 
cast in terms of the eigenvalues of this matrix. Eq. ([1]) 
can be written as x = — £x and, given that C is sym- 
metric, its solution in terms of normal modes is <pi(t) = 
(0)e~ Xit , where A 2 ; are the eigenvalues of £, see e.g. 
16j . The diffusion time scale r of the multiplex 



is controlled by the smallest non-zero eigenvalue of C. 
Specifically, r = I/A2. To get a physical insight on these 
eigenvalues as a function of the different diffusion coef- 
ficients within layers {D\ and D 2 ) and between layers 
(D x ), we propose to analyze the whole system using per- 
turbation theory. To simplify the notation, we choose 
the diffusion coefficients D\ = D 2 = 1 fixing then the 
relative time scale of the problem. 

Let us consider the decomposition C = Co + V, where 
Co is the block diagonal matrix corresponding to the 
Laplacians of every layer, with zeros in the off-diagonal 
blocks, and T> is formed by the rest of the elements. In 
matrix form it reads: 



C = Co + V = 




+ D X 












(3) 
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= D K Y,w^{xf -xf) + Y, D KL{x^ -xf) , (1) \\< ...A* and = \\ < A| < . . . A* , respectively, 



L=l 



where w^j denotes the weight matrix at layer K (wf^ = 
means that there is no link between nodes i and j in layer 
K) . This set of equations can be dimensionally lifted to 
a space of N x M dimensions. To have a more clear 
picture of our formalism wc will consider, without loss of 
generalization, the most simple case of two layers M = 2. 
First, wc define a column vector state of 2 AT elements, 



(x\ 



\ X \ 



- N 



) = (x^x 2 ) = x. Then Eq. (JT) can 



be written in matrix form, where the interaction matrix 
has a block structure that conforms an object we call 
supra-Laplacian C, with the same properties that any 
zero-sum rows Laplacian has: 



Let us start the discussion by considering D x = 0. 
Then, the eigenvalues of C are the set formed by the 
union of the eigenvalues corresponding to the Laplacians 
of each layer L\ and L 2 . The eigenvalues arc = A J < 
A 2 

while the eigenvalues of the supra-Laplacian matrix are 
= Ai = A 2 < A 3 < . . . < A 2 „, being A 3 = mm(A2, A|). 
It is interesting to note that to analyze the eigenvector 
space it is convenient to move to a new basis where the 
space corresponding to Ai = A 2 = is spanned by vec- 
tors (1 ■ • • 1 1 1 ■ • • 1) and (1 • ■ • 1| - 1 1) instead of the 

canonical (1 ■ ■ ■ 1|0 • • • 0) and (0 ■ ■ • 0|1 • • ■ 1). 

Now let us consider that the diffusion between layers 
is different from zero, D x ^ 0. In this case, the supra- 
Laplacian will have the trivial eigenvalue Ai = with cor- 
responding eigenvector (1---1|1---1), and a non-trivial 
eigenvalue A = 2D X that corresponds exactly to the 
eigenvector (1 • • • 1| — 1 • • • — 1), since 
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where L\ and L 2 are the respective Laplacians of each 
layer, and / is the identity matrix. Here we have replaced 
-D12 by D x to emphasize the role of the diffusion process 
among the same node at different layers. The Laplacian 
matrix of each layer K is just Lk — Sr — Wk, where 
Wr is the weights matrix at layer K, and Sk a diagonal 
matrix containing the strength of each node i at layer K, 
(Sk)h = sf — J2j w fj ■ Note that the diagonal block 
structure of the supra-Laplacian reflects the interaction 
within layers and the off-diagonal blocks the connectivity 
between layers. 



Note that this eigenvalue always exists, but it will be 
the smallest non-zero one only when D x is very small, as 
compared to Di and D 2 . 

Next, we focus our attention on the opposite limit, a 
very large diffusion coefficient between layers D x ^> 
1. Defining D x = 1/e, we can write 



C = D X 




D X C. (5) 



The spectrum of C is considered here a perturbation of 
that at e = 0. It is worth recalling that, for e = 0, the 
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spectrum corresponds to that of the coupling matrix 












(6) 



which consists of two eigenvalues (Ai =0 and A 2 = 2) 
both TV-degenerate and spanned by eigenvectors of the 
form (u|u) and (u| — u). i.e. vectors having identical or 
opposite values in the i th and (i + N) th components, re- 
spectively. Thus, in the limit D x — )• 00, the set of eigen- 
values of C will split in two groups, one showing a linear 
divergent behavior A « 2D X for the sub-space (u| — u), 
and another having a finite value A as the result of the 
undetermined limit (0 • 00) in Eq. ([5]) for the sub-space 
(u|u). 

Now, we use the common ansatz in perturbation the- 
ory and propose a perturbed solution in terms of both 
eigenvalues and eigenvectors: 



Xi 



- \(°) 



XT +e# ) +()(<?) 



(7) 



(i) 



where the super-indices within parentheses represent the 
order of the perturbation [13, EH • Given that a set of 
eigenvalues of C will diverge linearly as 2D Xl we concen- 
trate in proposing perturbations for the finite solutions. 
These correspond to consider the following perturbation 
of the eigenspectrum of C: 



A = + eX' , 



u 



u 



(8) 



Expanding to 0(e) the eigenvalue problem £v = Av we 
obtain: 

V K - u x ) + l 2 u y \U J 

Matching each of the components in Eq. (j9|) we get: 



Liu+ (ui - u' 2 )) = A'u , 
L 2 u+(u 2 -ui)) = A'u , 



(10) 



that, after adding and subtracting Eqs. (fT0|) . transform 
into: 



(L!+L 2 )u = 2A'u 

(Li - L 2 )u = 2(u; - u 2 )) 



(11) 



From the system of Eqs. (fTT|) it is revealed that u is an 
eigenvector of the network formed by the superposition 
of both layers' laplacians, and that the eigenvalue of C, 
at first order in the expansion, is 



A' 



2 



(12) 



being X s the eigenvalue of the superposition (L\ + L 2 ) 
corresponding to the eigenvector u. Moreover, given that 
the vector perturbation in Eq. (j8|) must be orthogonal 
(u|u) _L (u'jju^), we can also find the eigenvector of the 
superposition (L1 + L2) such that u 2 = — = — u', then 



u' = -(L 2 - Li)u. 



(13) 



Summarizing, the eigenvectors with finite (i.e. non di- 
vergent) eigenvalues of the supra-Laplacian C for a large 
value of the diffusion coefficient D x = 1/e between layers 
are 



eu 



with eigenvalue 



K 
2 



(14) 



being u and X s the eigenvectors and corresponding eigen- 
values of the superposition (L\ + L 2 ). 

The physical insight obtained is the following, for low 
values of the diffusion coefficient between layers, the dif- 
fusion time scale of the global system is controlled by 
the inverse of 2D X . This asymptotic result is valid until 
the order of D x is similar to those of D\ and Z? 2 . For 
large values of D x the eigenspectrum splits into a set of 
values that diverges as 2D X , and a set of finite values, 
associated to the superposition of the layers. The min- 
imal eigenvalue different from zero turns out to be half 
the eigenvalue corresponding to the superposition of both 
layers A s /2. 

A comparison between the diffusion time scale of the 
independent layers and the whole multiplex is possible 
using known bounds about the eigenvalues of the Lapla- 
cians [l9| | . The time scale associated to the multiplex 
for D x <C 1 is t = JIT" that means that the cross- 
diffusion between layers is the limiting value of the dif- 
fusion spreading. On the other hand, the time scale as- 
sociated to the multiplex for D x S> 1 is r w 2/A s . This 
latter case is far less trivial than the previous one. Using 
the bounds in [l9[ we deduce the following result: 



As 
2 



> 



A 



> min (A 2 , A 



(15) 



The above inequality implies that the diffusion in the 
multiplex will be faster than the diffusion in the slowest 
layer. Thus, as a consequence of the multiplex structure, 
at least one layer (the one with the larger diffusion time 
scale) has its diffusion speeded up. The emergence of a 
super-diffussion, i.e. the fact that the time scale of the 
multiplex is faster than that of both layers acting sepa- 
rately is, in general, not guaranteed and depends on the 
specific structures coupled together. Furthermore, the 
following inequality also holds (l9j : 



A s 2N 

— < 

2 ~ 2N - 



(16) 



being sf the strength of node i at layer K . 
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FIG. 2: (color online) Evolution of the eigenspectra of the toy 
model presented in Fig. [1] as a function of the coupling D x 
(top) for D\ = D 2 — 1, and comparison between the second 
smallest eigenvalues A2 of the different Laplacians (bottom) . 



Finally, it is worth noticing that although the previous 
analysis assumes that the networks within layers are con- 
nected, we have also analyzed the case in which this hy- 
pothesis is relaxed. Imagine for example two layers such 
as one layer has two disconnected components. In this 
situation, the results hold in the limit D x ^> 1, and in the 
limit D x <C 1 the lowest (different from zero) eigenvalue 
scales as aD x , with < a < 2 although the perturbed 
eigenvector is far more intricate. 

To illustrate our results, we have computed the evo- 
lution of the eigenvalues of the supra-Laplacian for the 
example represented in Fig. [TJ which corresponds to two 
random networks of N = 6 nodes. In Fig. [2] (top) we plot 
the eigenvalues as a function of the diffusion coefficient 
D x . We observe the splitting of the eigenvalues into two 
groups, divergent and finite values, as predicted. Fig. [2] 
(bottom) shows the theoretical estimates for A2 in the 
asymptotic limits D. x <C 1 and D x 1. Note that, 
except for the intermediate zone {D x ps 1), where the 



analysis does not hold, the agreement is excellent. In this 
panel we have represented, as indicated in the legend, the 
eigenvalues of each layer, the eigenvalue of the superposi- 
tion of both layers and the line corresponding to 2D X as 
well as the eigenvalue of the supra-Laplacian. The results 
undoubtedly confirm that both theoretical limits (small 
and large D x ) are correctly identified by the analytical 
derivations. Moreover, in this case we observe a super- 
diffusion process, being the time scale associated to the 
supra-Laplacian smaller than that of layer 1 and layer 2 
if they were considered independently, i.e. r < T\ < T<i. 

In conclusion, we have developed a formalism to un- 
veil the time scales of diffusive processes on multiplex 
networks. The approach has been specifically presented 
for a two-layer multiplex, in a particular set up in which 
nodes are preserved through layers. We obtained ana- 
lytical results in the two asymptotic limits of small and 
large diffusion coefficients between layers. The findings 
show that the multiplex structure is able to speed up the 
less diffusive of the layers. In principle, it could also give 
rise to a super-diffusion process thus enhancing the dif- 
fusion of both layers. This striking result appears when 
one considers that the diffusion between the layers of the 
multiplex is faster than that occurring within each of 
the layers. Thus, it paves the way to the analysis of 
super-diffusion processes in real multiplex scenarios such 
as multimodal transportation systems. On more general 
grounds, given the wide applicability of the properties 
of the Laplacian to address many dynamical properties 
of networked systems, our results constitute a first step 
towards a better understanding of linear and nonlinear 
processes on top of multiplex structures. 
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